Complicate dynamical properties of a discrete slow-fast predator-prey model with ratio-dependent functional response

Using a semidiscretization method, we derive in this paper a discrete slow-fast predator-prey system with ratio-dependent functional response. First of all, a detailed study for the local stability of fixed points of the system is obtained by invoking an important lemma. In addition, by utilizing the center manifold theorem and the bifurcation theory some sufficient conditions are obtained for the transcritical bifurcation and Neimark-Sacker bifurcation of this system to occur. Finally, with the use of Matlab software, numerical simulations are carried out to illustrate the corresponding theoretical results and reveal some new dynamics of the system. Our results clearly demonstrate that the system is very sensitive to its fast time scale parameter variable.


Complicate dynamical properties of a discrete slow-fast predator-prey model with ratio-dependent functional response Xianyi Li * & Jiange Dong
Using a semidiscretization method, we derive in this paper a discrete slow-fast predator-prey system with ratio-dependent functional response.First of all, a detailed study for the local stability of fixed points of the system is obtained by invoking an important lemma.In addition, by utilizing the center manifold theorem and the bifurcation theory some sufficient conditions are obtained for the transcritical bifurcation and Neimark-Sacker bifurcation of this system to occur.Finally, with the use of Matlab software, numerical simulations are carried out to illustrate the corresponding theoretical results and reveal some new dynamics of the system.Our results clearly demonstrate that the system is very sensitive to its fast time scale parameter variable.
In recent decades, human beings are suffering from some disasters of destroying the natural environment, such as pollution, species extinction, virus epidemics, etc.It is important to provide strategies to relieve the environmental pressure.Mathematical modeling can reveal the changing trend of the natural environment, therefore, more and more biologists, ecologists and mathematicians are committed to studying ecological balance using mathematical models.
The classical prey-predator model with prey-dependent functional response as proposed by Holling in 1965 1 is given by a system of coupled ordinary differential equations: where X and Y respectively denote the prey and predator densities at time T.Here both species are assumed to be distributed homogeneously within their habitats.The function R(X) represents prey's per capita growth rate.In this model, the prey-predator interaction is described by the prey-dependent function P(X) known as functional response, which quantifies the average amount of prey consumed by a single predator per unit of time.The predator is assumed to be specialist.The parameter e(0 < e < 1) known as the conversion efficiency determines the fraction of prey biomass that contributes to the predator's growth.The function M(Y) represents per capita death rate of predators in the absence of prey.The important assumptions of the classical model are that predator encounters prey at random and the trophic function depends on prey abundance only.However, in the late 1989, Adriti and Ginzburg 2 challenged the classical theory by showing the importance of predator interference whenever the prey abundance is low.The authors argued that "the trophic function must be considered on the slow time scale of population dynamics at which the models operate-not on the fast behavioral time scale" 2 .It is therefore reasonable to assume that the functional response depends on the ratio of prey to predator abundance rather than just on the prey abundance when the available prey density is low.That is, in order to reflect the predator interference, the per capita functional response should be a function of X/Y rather than X.Based upon this idea, the Michaelis-Menten-Holling type functional response was introduced, also known as www.nature.com/scientificreports/ratio-dependent functional response 2 .In this paper we consider the ratio-dependent prey-predator system with the logistic growth in prey as a baseline model as follows where r is the linear growth rate of the prey, k is the environment carrying capacity to prey, b is the maximum per capita growth rate of the predator (note that the corresponding term approaches its limiting value bY when X becomes very large), m is the predator mortality and g is the relative saturation factor between the two species.
By rescaling the variables we obtain the following dimensionless system where the new parameters ǫ = b r , α = a gr and δ = m b are dimensionless.Generally, one assumes that the prey population grows faster than the predator population (which, in fact, is often the case in nature).So, we have b < r , implying 0 < ǫ < 1.
Note that the dimensionless time τ in the system (3) is the slow time.With the transformation t = τ ǫ , the equivalent system in fast time scale is It is not easy to solve a complicate differential equation (system) without computer.So, one tries to use discretization method to derive and study the discrete model of a complicate differential equation (system) so that one can better understand the properties of corresponding continuous system.How to discretize a complicate differential equation (system)?Many discretization methods, such as the forward Euler method, the backward Euler method, the semidiscretization method, and so on, can be utilized.The discrete version of the system (4) has not been investigated yet.We here use the semidiscretization method to study its discrete model.The advantage for this kind of discrete method is for one not to require to consider the step size.Relatively speaking, this kind of method can reduce the number of parameters so that the system studied is easily investigated.
For this, suppose [t] to denote the greatest integer not exceeding t.Consider the average change rate of the system (4) at integer number points It is easy to see that the system (5) has piecewise constant arguments, and that the solution (x(t), y(t)) of the system (5) for t ∈ [0, +∞) possesses the following characteristics: 1. on the interval [0, +∞) , x(t) and y(t) are continuous; 2. when t ∈ [0, +∞) except for the points t ∈ {0, 1, 2, 3, . ..} , dx(t) dt and dy(t) dt exist everywhere.
Letting t → (n + 1) − in (6) produces where the parameters a > 0, δ > 0, 0 < ǫ < 1 have the same meanings as in the system (4).We mainly consider in this paper the dynamical properties of the system (7).The rest of this paper is organized as follows: In section "Existence and stability of fixed points", we investigate the existence and stability of fixed points of the system (7).In section "Bifurcation analysis", we derive the sufficient conditions for the transcritical bifurcation and the Neimark-Sacker bifurcation of the system (7) to (2) (5) occur.In section "Numerical simulation", numerical simulations are performed to illustrate the theoretical results derived and reveal some new dynamical properties of the system.
Before we analyze the fixed points of the system (7), we recall the following lemma 3 .
, where B and C are two real constants.Suppose 1 and 2 are two roots of F( ) = 0 .Then the following statements hold.
1 and 2 are a pair of conjugate complex roots and,

Existence and stability of fixed points
In this section, we first consider the existence of fixed points and then analyze the local stability of each fixed point of the system (7).
The fixed points of the system (7) satisfy Considering the biological meanings of the system (7), one only takes into account nonnegative fixed points.Thereout, one notices that the system (7) has and only has three nonnegative fixed points The Jacobian matrix of the system (7) at any fixed point E(x, y) takes the following form The characteristic polynomial of Jacobian matrix J(E) reads where For the stability of fixed points E 0 , E 1 and E 2 , we can easily get the following Theorems 1, 2 and 3 respectively.

Theorem 2
The following statements about the fixed point E 1 = (1, 0) of the system (7) are true.
The proofs for Theorems 1 and 2 are easy and omitted here.
) is a positive fixed point of the system (7).
and δ 0 be the unique positive root of the function f (δ) = δ 3 − δ 2 + δ − a−1 a for a > 1 and δ ∈ ( a−1 a , a−1 a ) .Then the following statements are true about the positive fixed point E 2 .

Bifurcation analysis
In this section, we use the center manifold theorem and bifurcation theory to analyze the local bifurcation problems of the system (7) at the fixed points E 1 and E 2 , rspectively.For related work, refer to 6-12 .
Theorem 4 Suppose the parameters (a, δ, ǫ) ∈ S E 1 .Let δ 0 = r a+1 , then the system (7) undergoes a transcritical bifurcation at the fixed point E 1 when the parameter δ varies in a small neighborhood of the critical value δ 0 .
Proof In order to show the detailed process, we proceed according to the following steps.
Let u n = x n − 1, v n = y n − 0 , which transforms the fixed point E 1 = (1, 0) to the origin O(0, 0), and the system (7) to Giving a small perturbation δ * of the parameter δ , i.e., δ * = δ − δ 0 , with 0 < |δ * | ≪ 1 , the system ( 8) is perturbed into Letting δ * n+1 = δ * n = δ * , the system (9) can be written as Taylor expanding of the system (10) at (u n , v n , δ * n ) = (0, 0, 0) takes the form where Taking the transformation   4 , or refer to 5 , all the conditions for the occurrence of a transcritical bifurcation are established, hence, there is an occurrence of a transcritical bifurcation for the system (7)  at the fixed point E 1 .The proof is over.

Bifurcation at E 2 -Neimark-Sacker bifurcation
First from the proof process of Theorem 3, we know that F(1) > 0 and F(−1) > 0 .Namely, 1 and −1 are not eigenvalues of the system (7) at the fixed point E 2 .So, it is impossible for the system (7) at the fixed point E 2 to undergo a transcritical bifurcation, or a fold bifurcation, or a pitchfork bifurcation or a flip bifurcation.
3 ).www.nature.com/scientificreports/But, from Theorem 3 one also sees that, when a > 1 , 2 , E 2 is non-hyper- bolic.Moreover, when the parameter ǫ crosses the critical value ǫ 0 , the dimensional numbers of the stable mani- fold and the unstable manifold of the system (7) at the fixed point E 2 have an essential change.So, a bifurcation will occur at this time.Indeed, we derive that the system (7) at the fixed point E 2 will undergo a Neimark-Sacker bifurcation in the space of parameters where δ 0 is the unique positive root of the function f (δ) = δ 3 − δ 2 + δ − a−1 a for a > 1 and δ ∈ ( a−1 a , a−1 a ).In order to show the process clearly, we carry out the following steps.Take the change of variables u n = x n − x 0 , v n = y n − y 0 , transforming the fixed point E 2 = (x 0 , y 0 ) to the origin O(0, 0) and the system (7) into Give a small perturbation ǫ * of the parameter ǫ around the critical value ǫ 0 , i.e., ǫ * = ǫ − ǫ 0 , then the perturba- tion of the system (13) can be regarded as follows The corresponding characteristic equation of the linearized equation of the system (14) at the equilibrium point (0,0) can be expressed as where and Noticing p(ǫ * ) ǫ * =0 = a(1 − δ 2 ) + 1 a(1−δ) − δ and q(ǫ * ) ǫ * =0 = 1 , one finds that p 2 (0) − 4q(0) < 0 always holds.In fact, it is easy to see It follows from a > 0 and δ ∈ (0, . This is verified by (a, δ, ǫ) ∈ S E + .So, when 0 < |δ * | ≪ 1 , the two roots of F( ) = 0 take as where i is an imaginary unit, namely, i 2 = −1 ; moreover which implies The occurrence of Neimark-Sacker bifurcation requires the following two conditions to be satisfied The transversal condition H 1 obviously holds.For the sake of convenience in discussing the nondegenerate condition ( H 2 ), let Then 1,2 (0) = α 1 ± α 2 i.It is easy to derive m 1,2 (0) = 1 for all m = 1, 2, 3, 4 .Hence, (H.1) and (H.2) hold.Accord- ing to [4, page 517-522], all the conditions for a Neimark-Sacker bifurcation to occur are satisfied.
In order to derive the normal form of the system (14), we expand the system (14) in power series to the third-order term at the origin to obtain where Make a change of variables (u, v) T = T(X, Y ) T with matrix T = a 01 0 α 1 − a 10 − α 2 , then the system (15) is transformed as where  www.nature.com/scientificreports/Our results also clearly demonstrate that the system ( 7) is very sensitive to its fast time scale parameter variable ǫ , namely, to appropriately adjust the value of fast time scale parameter variable ǫ may alter the stability of the system (7).So, the results in this paper also provide a way for how to control the stability of the system (7).